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O ; Abstract 

o 

"^ \ A simple model of epitaxial growth proposed by Wolf and Villain is inves- 

Y* , tigated using extensive computer simulations. We find an unexpectedly com- 

^ , plex crossover behavior of the original model in both 1+1 and 2+1 dimensions. 

O 

A crossover from the effective growth exponent /3eff~0.37 to /3cff~0.33 is ob- 
k^ , served in 1+1 dimensions, whereas additional crossovers, which we believe are 

5^ I to the scaling behavior of an Edwards- Wilkinson type, are observed in both 

1+1 and 2+1 dimensions. Anomalous scaling due to power-law growth of 
the average step height is found in 1+1 D, and also at short time and length 
scales in 2+1 D. The roughness exponents C,^^ obtained from the height-height 
correlation functions in 1+1 D (~3/4) and 2+1 D (~2/3) cannot be simul- 
taneously explained by any of the continuum equations proposed so far to 

describe epitaxial growth. 
68.55.-a,68.35.Fx,61.50.Cj,64.60.Ht 



Typeset using REVTgX 



Simple ( "toy" ) models of surface growth have become very popular recently, mainly as 
examples of externally driven non-equilibrium systems that exhibit simple (yet nontrivial) 
scaling behavior.EJ The focus is on the surface width (roughness) w defined asw = {\J h^ — h ) 
(the overbars denote spatial averages, () a statistical average), i.e., the variance of the surface 
height profile h{r,t). The surface width scales with the time t and the linear dimension of 
the substrate L as w(t,L) oc L''f(t/L^), where the scaling function /(x) has the properties 
/(x) = const, X 3> 1 and /(x) oc x'^,x <^ 1, /5 = (/z. Thus, w grows according to a power 
law, w oc t^, until a steady state characterized by a constant value of the width is reached 
after a time tgat proportional to L^. The value of the saturated width Wsat varies with the 
system size according to Wsat oc L"^. The exponents ( and (3 (or ( and z) characterize the 
scaling behavior of the roughness for a particular model and determine its universality class 
in analogy with theory of critical phenomena. 

Alternatively, one can study the surface roughness using the height-height correla- 
tion function G{v,t) = ([/i(x + r, t) — /;,(x, t)]^) which obeys the scaling relational G{r,t) oc 
r"^^" g{r /t^/^''), where the scaling function g{x) is constant for x <^ 1 and g{x) oc x"^''" for 
x^^l (equivalently, the structure factor S can be used, see e.g. Ref. |^). In most of growth 
models the exponents obtained using the two different methods are equal.EJ 

Scaling behavior of the surface roughness can be investigated using continuum, stochastic 
differential equations which in the case of conserving models with surface diffusion have the 
form 

^^ = -V.j(r,t)+r^(r,t) , (1) 

where ?7(r, t) is a zero mean, random noise term in the incoming flux, and the current j(r, t) is 
a function of the derivatives of /i(r, t). In a number of recent theoretical studies,ErEJ models 
in which surface diffusion is the dominant physical mechanism of the surface smoothing 
were studied. The scaling relation 2( = z ~ d' (where d' is the substrate dimension) holds 
for these models.Q The most often studied cases were j oc — V/i [Edwards- Wilkinson (EW) 
inodelEj], j oc VV^/i (the linear diffusion model&EJ), and j oc V(V/i)^ and j oc (V/i)^ (the 



nonlinear diffusion inodelsQlj'tll which we will denote I and II, respectively). The predicted 
values of exponents arJ'li (3^^ = (3 - d)/A, C^^ = (3 - d)/2, (3^^''' = (5 - d)/8, ('*" = 

(5 - rf)/2, iS-^onlin-I^^^ _ ^y^^ ^ ^)^ ^nonlin-I^^^ _ ^y^^ ^^^ pnonlin- 1 1 ^ (^^ _ ^)/2(3 + d) , 

^nonhn-ii _ ('5 _ c?)/4 where d = d' + 1. It should be noted that these values are usually 
based on renormalization group (RG) calculations within the one-loop approximation (e.g., 
Refs. |6| Jl^) or on Flory-type approximation (Ref. ^. 

Alternative approach is to employ a powerful computer and study discrete models with 
microscopic rules reflecting physically important surface processes.EJ From a wide variety of 
discrete models we focus our attention on the one proposed by Wolf and VillainQ (WV) which 
is based on the solid-on-solid model of crystal growths and has been designed with growth 
by molecular-beam epitaxy (MBE) in mind. The growth rules of the model (see below) are 
supposed to mimic surface diffusion, the principal mechanism of surface smoothing during 
MBE. Original simulations of the modela in 1+1 D yielded exponents /3eff = 0.365±0.015 and 
Ceff = l-4±0.1 (thus Zeflf = 3.8±0.5) in agreement with the theoretical prediction of the linear 
model. However, a subsequent numerical workQ has shown that in 2+1 D the values of the 
exponents are /3eff = 0.206+0.02 and Cefr = 0.66+0.03 (thus ^efr = 3.2 + 0.5) which correspond 
to the prediction of the nonlinear model I. The puzzling difference between the behavior of 
the model in 1+1 and 2+1 D has been confirmed in another numerical study.0 The possible 
explanation could be a slow crossover from exponent of the linear model to the exponent 
of the nonlinear model I observable in 1+1 D only for large system sizes and long times, 
similar to a crossover observed in a full-diffusion solid-on-solid model.B To complicate the 
situation even further, it has been proposed recently by Krug, Plischke, and SiegertEj that 
the WV model should cross over to EW behavior for more than ^ 10^ (1+1 D) or ?^ 2 x 10'^ 
(2+1 D) deposited layers. This result was based on the study of the inclination-dependent 
diffusion current which is supposed to generate the EW term (V^/i) in continuum differential 
equations. t3'tll The EW term is more relevant (in the RG sense) than all allowed nonlinear 
termsti and governs the asymptotical behavior of the model. The long time needed to 
observe the asymptotic regime may be explained as due to a very small coefficient in front 



of the EW term. Finally, it has been found very recently that the WV model in 1+1 D 
does not fulfill standard scaling and that different values of the exponents are obtained from 
behavior of the surface width and the correlation function or the structure factor.liil The 
exponents obtained in 1+1 D from the behavior of the structure factor were Cgg = 0.75+0.05 
and 2;gfj = 2.4 + 0.1 (thus /?gg~0.31) very close to the nonlinear model II predictions.EJ The 
authors also suggested that this anomalous behavior is present in 2+1 D as well. 

The need for better simulation data apparent from the above summary has motivated 
the present work. We performed large-scale simulations of the WV model in 1+1 and 
2+1 D. In 1+1 D, we used system sizes as large as L = 40 000 sites and deposited up to 
2^^ ^ 1.3 X 10^ layers whereas in 2+1 D we used lattice sizes of up to 1000 x 1000 depositing 
up to 2^^ ^ 1.3 X 10^ layers. Our results for exponents obtained from the surface width 
show (i) a crossover from Pcs ~ 0.37 (/5''") to /?eff ~ 0.33 (^"o"'«"-^) in 1+1 D, (ii) crossovers 
from /?eflF ~ i^nonhn-i ^^ ^]^g scaliug bchavior which we believe is of the EW model in both 
1 + 1 and 2+1 D. Exponents calculated from the correlation function are in agreement with 
those of the nonlinear model II in 1+1 D [Q^ k, 0.75) and the nonlinear model I in 2+1 D 
(Cff ^ 0.65). 

The microscopic rules of the basic model are the same as in Ref. ^. In every time step, a 
particle is added at a randomly chosen lattice site and then relaxes toward a nearest-neighbor 
site which offers the highest coordination (the number of nearest neighbors) where it sticks 
for the rest of the simulation. If the number of nearest neighbors cannot be increased, 
particularly in the case of tie (one or more neighboring sites have the same coordination as 
the original site) the particle stays at the initial position. The generalization to 2+1 D is 
straightforward. 

In 1+1 D, we carried out simulations for lattice sizes L=150, 300, 600, 800, 1000, 2000, 
and 40 000 depositing from 2^^ to 2^^ layers (except for L = 40 000 where only 2^^ layers 
were deposited) , see Fig. |I]. (Note that the curves for lattice sizes L > 800 have been offset 
in Fig. |I| to avoid confusion due to overlapping data points; statistical errors of the data 
points are not larger than the symbols size.) The crossover of the exponent (3ea from the 



value 0.372 ± 0.007 corresponding to the linear model (/3'*" = |) to ~ 0.336 ± 0.009 in very- 
good agreement with the nonlinear model I prediction (^^"o"'«"-'f = 1) can be best observed 
for the largest simulated lattice, L = 40 000. It is apparent from Fig. |I] that this crossover 
takes place after approximately 2^^ ^ 3 x 10'^ layers have been deposited. The values of Wgat 
we have obtained (48 ± 1, 91 ± 10, and 209 ± 20 for L = 150, 300, and 600, respectively), give 
Ceff = 1-05±0.1, in very good agreement with the nonlinear model I prediction, (^^"'™-^ = i_ 

However, a more surprising observation can be made after close inspection of the data in 
Fig. |l[ All the curves (with the possible exception of the one for L = 40 000 where fewer layers 
were deposited) show another decrease in the value of the exponent P^s after the deposition 
of 2^^ ^2 X 10® layers. The curves for L>800 seem to follow the slope corresponding to the 
value /? = i of the EW model in 1+1 D (see Fig. p. Because the time tsat{L) needed for the 
saturation depends on the system size as L^, we obtain tsat(-^ = 2000) > 2^'^ based on our 
results for tsat for small lattice sizes (even if we use 2;"-°"'*""^^ = 5/2). Moreover, the position 
of the crossover seems not to depend on the lattice size. Hence we can exclude that this is 
the effect of the saturation. 

Similar behavior is revealed by an inspection of 2+1 D results in Fig. ^ The data given 
for the lattice sizes L = 500 and L = 1000, correspond at first (between ^ 2^ and 2^° layers 
deposited) to the value of PeS = 0.22 ± 0.05, very close to (but definitely higher than) the 
nonlinear model I prediction |. After ^ 10^ layers are deposited, a crossover takes place. 
Again it cannot be due to the vicinity of the saturation regime because the estimate of tgat 
in 2 + 1 Z) for L = 500 is larger than 2^^ layers based on our previous resultsQ and the position 
of the crossover does not shift with the lattice size. A logarithmic increase of the roughness 
(/3^^ = 0) is expected for the EW model in 2+1 D. Double-logarithmic and semilogarithmic 
plots of the last seven data points (averaged over the both curves) are shown in the inset 
of Fig. ^. The statistics of the data are unsufficient to draw a definitive conclusion, but it 
seems that the semilogarithmic plot follows a straight line whereas the log-log one is curved. 
Thus, we believe our data suggest that in both 1+1 and 2+1 D we observe crossovers to the 
scaling behavior of the EW model predicted by Krug et alx^ and it should be noted that 



even the positions of these crossovers agree well with this prediction. 

Following the paper by Schroeder et a/.lHl we also studied the behavior of G{r, t) . The time 
dependence of G(r, t) in 1+1 and 2+1 D is shown in Fig. ^. The exponent (^^g = 0.75+0.03 we 
obtained in 1+1 D using the initial slope of the correlation function for t = 2^^ monolayers and 
L = 40 000 (Fig. 1^ (a)) differs substantially from the value Coff~l obtained from the surface 
width behavior at later times (see above) and corresponds to (^"o»^'«"--f-'' = 3/4. However, the 
value Ceff = 0-65±0.03 obtained in 2+1 D (Fig. |] (b)) for L = 500 and t = 2^'' monolayers is 
definitely different from 2+1 D value of ("o'^'*"-^^ = 1/2 and is much closer to («°"'i"-^ = 2/3 
which was also obtained from the surface width behavior.Q The values of the exponents Q^ 
obtained from G{r^ t) are thus to the best of our knowledge inconsistent with any simple 
model of epitaxial growth proposed so far. 

In agreement with Ref . |Tl|, we see the initial power-law increase of the average step height 
G{l,t)oct^ with A = 0.38 + 0.01 in 1+1 D, and A = 0.095 + 0.003 in 2+1 D, see the insets in 
Fig. ^ In Ref. |Tl|, it is shown how this can explain effective exponents corresponding to the 
linear equation and the discrepancy between Qg and Ces i^ 1+1 D provided the nonlinear 
model II is used. Notice also that the crossover from the exponent /?eflf~ 0.375 to the value 
/3eflf ~0.33 coincides with the beginning of G(l, t) saturation. In 2+1 D, if we use the values 
of the nonlinear model I and A ~ 0.095, we obtain (es ~ 0.84 and (3es ~ 0.23. The value of 
/Seff is in agreement with the results obtained here for short times t < 2^° [see above; notice 
again that this is the region where we observe the power-law growth of G(l, t)] and also the 
value of Ceff is not inconsistent with the initial slope of the plot w vs. L in Fig. 2 of Ref. |^. 



Hence, the theory of Ref. |TT] agrees with our results provided we use different "underlying" 
models in 1+1 and 2+1 D. The anomalous behavior due to power-law growth of G{l,t) is 
much weaker in 2+1 D but can be observed at short time and length scales. 

We also tried to find changes in the correlation function behavior that should take place 
following the crossovers which we suppose are to the EW behavior (see the thick lines in 
Fig. I). In 1+1 D [Fig. | (a)], a decrease to the value of Ceff~2/3 is observed for L = 2000 
and 2^^ monolayers deposited, but this is still significantly higher than the EW model value 



(■^'^ = 1/2. In 2+1 D [Fig. ^ (b)], the behavior of the correlation function after the crossover 
does not change appreciably. A probable explanation is that the range of length scales (larger 
than the relevant crossover length) where G{r, t) is influenced by the EW behavior is too 
short to be resolved in our data. However, we have found that the structure factor calculated 
for L = 2000 and 2^^ monolayers does at large wavelengths follow the slope —2 expected for 
the EW behavior. Statistics of these data are not sufficient to provide unambiguous proof 
and more work is needed. 

In conclusion, we have studied kinetic roughening in 1+1 and 2+1 dimensions of a simple 
model of epitaxial growth proposed by Wolf and Villain.Q From the study of the surface width 
in 1 + 1 D, we have found a crossover from the scaling behavior of a linear differential equation 
proposed to describe MBE growth (/3eff ~0.37) to that corresponding to a nonlinear model I 
equation (/?eff ~0.33). According to Ref. |TT| and our results, unusual behavior of the model 
in 1+1 D is due to a power-law growth of the average step height, G{1, t) oct^ where A~0.38 
in 1+1 D and ~ 0.095 in 2+1 D. As a consequence, the exponent Qf^ obtained in 1+1 D 
from the height-height correlation function differs from the value obtained from the scaling 
of the surface width and does not change during the crossover. Its value in 1+1 D (Qq^S/A) 
corresponds to another nonlinear model II proposed by Lai and Das Sarma.1^ However, the 
value of the exponent Qg ~ 2/3 in 2+1 D is consistent with the nonlinear model I. One 
possibility of how to explain the discrepancy between the behavior in 1+1 and 2+1 D in 
terms of a continuum equation is to suppose that both relevant nonlinearities [V^(V/i)^ and 
V(V/i)^] are present but with different coefficients (and thus different importance) in 1+1 
and 2+1 D. 

In both 1+1 and 2+1 D, we see additional crossovers and our results seem to indicate 
these are to the scaling behavior of the EW model. The conclusion that the asymptotic 
behavior of the "ideal" MBE growth model is of EW type is in agreement with the recent 
suggestion by Krug et a/.t3. It would be interesting to see this crossover also in other 
discrete models, in particular in the full-diffusion model of Ref. ^ which was very successful 
in describing the initial stages of MBE growth. In some variants of the models with surface 



diffusion, the final crossover to the EW behavior may be absent and instead the instabihty 
may develop. We expect this kind of behavior in models with barriers to interlayer transport. 
We thank Professor D. Wolf and Professor D.D. Vvedensky for stimulating discussions 
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FIGURES 
FIG. 1. Surface width vs. time for the WV model in 1+1 D (L = 150(A), 300(x), 600 (d), 

800(o),1000(»),2000(o), and 40 000 (— )). Notice that data for larger lattice sizes (L > 800) were 

offset to avoid overlapping of data points. 

FIG. 2. Surface width vs. time for the WV model in 2+1 D (L = 500(A) and lOOO(o)). Data 
for L = 1000 were offset to avoid overlapping of data points. The inset shows the last seven data 
points for L = 1000 in semilogarithmic (•) and double logarithmic (o) scales (see text). 

FIG. 3. (a) The height -height correlation function G{r,t) in 1+1 D for L = 40 000 and the times 
t = 2i(o),24(.),27(o),2i0(x),2i3 (n), 2i6(*), 219(A), 222 (+), ^^^ ^^ L = 2000 and ^ = 2^7 (_); (b) 

The height -height correlation function G{r,t) in 2+1 D for L = 500 and the times t = 22(*),2^ (d), 
26(x),2^(A),2i°(»),2i2(o)^2i^(o), and 21^^ (— ). The insets show the power-law increase of G(l,t). 
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